Limits to squeezing in the degenerate OPO 
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We develop a systematic theory of quantum fluctuations in the driven parametric oscillator 
(OPO), including the region near threshold. This allows us to treat the limits imposed by nonlin- 
earities to quantum squeezing and noise reduction, in this non-equilibrium quantum phase-transition. 
In particular, we compute the squeezing spectrum near threshold, and calculate the optimum value. 
We find that the optimal noise reduction occurs at different driving fields, depending on the ratio 
of damping rates. The largest spectral noise reductions are predicted to occur with a very high-Q 
second-harmonic cavity. Our analytic results agree well with stochastic numerical simulations. We 
also compare the results obtained in the positive-P representation, as a fully quantum mechanical 
calculation, with the truncated Wigner phase space equation, also known as semiclassical theory. 
PACS numbers: 03.65.Bz 



I. INTRODUCTION 



Optical parametric oscillators are one of the most interesting and well characterized devices in nonlinear quantum 
£NJ , optics. Novel discoveries made with them include demonstrations of large amounts of squeezing |l[ , significant quantum 
intensity correlations together with a quadrature correlation measurement that provided the first experimental 
demonstration H of the original EPR paradox. Practical applications include their use as highly efficient and tunable 
frequency converters. In the present paper, we focus on the optimum below-threshold squeezing results, which 
determine the limits to squeezing obtained near the critical point, where nonlinear corrections start to dominate. In 
a companion paper Q| , the related question of critical fluctuations at threshold is treated. 
iy) ' The theory of quantum squeezing in the linear parametric oscillator is now well-developed E3|. Excellent 
agreement between theory and experiment is obtained 0, in the region below threshold. However, the usual 
theory is linearized, and therefore cannot be used in the near-threshold region where the squeezing is largest. The 
drawback with linearized theories is that they predict that zero quantum noise levels are achievable at threshold. 
■ This is clearly unrealistic, since (by the Heisenberg uncertainty principle) it necessarily requires an infinite energy 
fi i| in the conjugate mode. More significantly, this would imply an infinite amount of phase information - which is also 
, 1 * . impossible, since the coherent pump which drives the parametric oscillator can only supply a finite quantity of phase 
£h ' information. 

While present experiments are limited by technical noise from approaching the critical point too closely, it is 
reasonable to expect that progress in integrated optics will lead to more stable, highly miniaturized devices which 
could well operate at the quantum limit, even near threshold. Accordingly, there have been a number of investigations 
as to the ultimate limits to the squeezing spectrum of a parametric amplifier/oscillator. This has often involved using 
Keldysh diagrams or Wyld-Keldysh techniques |l(J]- Jl2| to extend the linear theory jl3|]- JL5J, using a many-body 
theory analog of Feynman diagrams. 

The two-mode Hilbert space involved in these problems typically has a minimum dimension > 10 6 , even with only 
N = 10 3 photons , and therefore would be difficult to solve using other methods that involve number-state expansions 
- either using a direct solution of the master equation, or stochastic wave-function methods. The Hamiltonian matrix 
would have 10 12 coefficients, unless simplified, with a density matrix of similar size. More typical experimental photon 
numbers have at least N = 10 6 , with a corresponding density matrix dimension of 10 24 - which appears completely 
inaccessible with number state techniques. Another drawback of number-state techniques is that they usually do not 
permit analytic approximations, which can give more physical insight. 

We therefore treat these questions using the coherent-state positive P-representation p6| , combined with an expan- 
sion technique valid below the critical point. Results are also verified by the use of direct numerical stochastic equation 
simulations. We find an N~ 2 / 3 scaling law for the optimal squeezing predicted by Plimak and Walls H] is obtained 
here as well, but with a different spectrum, owing to the use of more systematic expansion techniques that result 
from using the positive P-representation method. Our analytic results for optimal squeezing, which occurs below 
the critical threshold, give excellent agreement with accurate numerical simulations for the same parameter values. 
However, even larger noise reductions are predicted to occur simply by reducing the losses of the second-harmonic, 
in which case the N~ 2 ^ 3 scaling law no longer holds. In a companion paper we consider the related problem of the 
critical region, where the narrow-band squeezing is less than optimal due to the effects of critical fluctuations. 
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We also compare the above results with a semiclassical approach, that is, a truncated Wigner phase-space equation. 
This equation corresponds to a classical theory with added vacuum fluctuations. A comparison between the positive 
P-representation (fully quantum mechanical) and semiclassical theories permits us to see how far one can go and what 
is the limitation of this extended classical point of view. We find that the nonlinear corrections in the semiclassical 
theory are in strong disagreement with the full quantum theory far below threshold, but agree near threshold. This 
tells us that the semiclassical theory works surprisingly well in the threshold region, indicating that the large quantum 
fluctuations near threshold have a rather classical character. 



II. HAMILTONIAN AND MASTER EQUATION 



The model considered here is the degenerate parametric oscillator. The system of interest is an idealized interferom- 
eter, which is resonant at two frequencies, lo\ and to 2 « 2u)%. It is externally driven at the larger of the two frequencies. 
Both frequencies are damped due to cavity losses. Down conversion of the pump photons to resonant sub-harmonic 
mode photons occurs due to a nonlinearity present inside the cavity. The Heisenberg picture Hamiltonian that 
describes this open system || is 



H = H. 



sys 



J2 fc^ft + otr 



3=1,2 



3 3 



H, 



(2.1) 



where the intra-cavity or system Hamiltonian is given by: 
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Here £ represents the external driving field at frequency ui 2 . The term Hr describes the free evolution of the 
extra-cavity modes that are the loss-reservoirs of the cavity. The term \ 1S the coupling parameter due to a \^ 
nonlinear medium internal to the cavity, and , Tj are reservoir operators that create and destroy photons in the 
loss-reservoir coupled to the internal mode of frequency u>j . 

Next, we wish to consider an interaction picture, obtained with the definition that 



Hq = hujja-a,j 



(2.3) 



In other words, the operators will evolve according to the relevant mode frequency, while the states evolve according 
to the rest of the system Hamiltonian. The interaction Hamiltonian used here then reduces to the the standard one 
M for a non-degenerate, single-mode parametric amplifier or oscillator: 
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(2.4) 



Here ai,a 2 are now time- independent operators representing the fundamental and second- harmonic modes respec- 
tively. For simplicity, we have chosen the field mode-functions so that £, x are real. 

Using standard techniques JTt]] to eliminate the heat bath, we obtain the following master equation for the reduced 
density operator of the system in the interaction picture: 
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(2.5) 



where ji are the internal mode amplitude damping rates and we assume that n,; << 1 , where ni are the mean numbers 
of thermal photons in the input reservoir modes. Using reservoir theory, it is possible to identify the coherent driving 
field with a corresponding input photon flux from an external coherent laser, with I 2 = \£\ 2 /2j2 n in photons / s , where 
7™ is the input coupler decay rate. For optimum performance, we will assume that 72 = 72™ , and similarly for the 
fundamental mode - which will be assumed to only decay through its output coupling mirror. If these conditions are 
not satisfied, then the coupling efficiency and maximum squeezing arc reduced. 

At this point, we note that in the classical limit, the system has the well-known classical equations of intra-cavity 
parametric oscillation, where we define: on — (a^), and hence obtain, in the interaction picture: 
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(2.6) 



These equations are valid in the limit of large photon number. They are obtained by the use of a classical decorrelation 
in which all operator products are assumed to factorize, so that (o|ftj) — {a\)(aj), and (ffliflj) — (a,i)(aj). The solution 
of these equations is immediate classically, and has the property that there is a phase transition at the critical driving 
field of: £ = £ c = 7172 /x ■ F° r driving fields below this value, one has: 



a x = 0, 
a 2 = £/7a ■ 

For fields above this value, the signal field ct\ is bistable, with: 



(2.7) 
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71 

«2 = ■ 
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(2.8) 



The intra-cavity photon number at the critical point is N c — jf/x 2 — £ 2 M2 ■ Classically, there are only second- 
harmonic photons present at this driving field, and the input photon flux is I c = £^/2-f 2 = 7i72/4x 2 - However, a 
squeezed field - with finite intensity - is actually emitted as well. This is not taken into account in the classical theory. 



III. OPERATOR REPRESENTATIONS 



In order to treat the full quantum evolution, we now turn to the methods of operator representation theory. These 
techniques can be used to transform the density matrix equations of motion to c-number Fokker-Planck or stochastic 
equations. 



A. The positive P-representation 

In the positive P-representation, the density matrix is expanded in terms of multi-mode coherent state vectors | a?): 

p= / P(a, a^) -j—^r* ' d ad . (3.1) 
J ((a^) \a) 

Following standard procedures, the assumption of vanishing boundary terms allows the master equation to be 
re-written as a Fokker-Planck equation in P (c?73 + ) , and hence as a stochastic equation fl§|| with real noise. The 
assumption of vanishing boundary terms is critical to this procedure, and we note here that this is generally valid 
when the ratio of nonlinearity to damping is small |ist| , (i.e. |x/7fe| ^ !)• The stochastic procedure is best regarded as 
being generally an asymptotic procedure, valid for small |x/7fc| " m which case the boundary terms are exponentially 
suppressed. We check this assumption numerically here as well, and point out that the required ratio of nonlinearity to 
damping is extremely well-satisfied in current experiments, where the ratio is typically 1CP 6 or less. Further analysis 
of thi s pro blem has been given elsewhere . Given this assumption, the following stochastic equations are obtained 
from ( |2.5| ) and ([O]), for any driving field £, that is, either below or above threshold: 



doti = 
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The stochastic correlations are given by: 

(dw k (t)) = , 

(dw k (t)dwi(t)) =6 kl dt . (3.3) 

This means that dwk(t) represent two real Gaussian and uncorrected stochastic processes, and the amplitude of 
the stochastic fluctuations that act on the signal mode are dependent on the pump field dynamics. Our derivation is 
formally based on the Ito stochastic calculus. However, in this case, either Ito or Stratonovic stochastic calculus gives 
identical results [Q . 



B. The semiclassical theory 



We can also write a c-number phase space equation using an approximate form of the Wigner representation jS| , which 
is equivalent to stochastic electrodynamics. The characteristic function of the Wigner representation is written as 



Xw(z) = Tripe 



iz* a) -\-iza 



Tr ( pe 



a iz* a' p iza p — \ z\ 2 /2 



and the Wigner distribution can be written as Fourier transform of the characteristic function 



W(a) 



1 



d 2 zxw( z )e 



(3.4) 
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In the Wigner representation, the phase space equation that corresponds to the master equation (2.5) is 
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(3.6) 



If we truncate the third derivative of the phase space equation we get a genuine Fokker-Planck type equation with 
positive definite diffusion constant. This can be mapped into the following Ito stochastic differential coupled equations 
(for simplicity we let fii — 0, as before). 

dai = [—710:1 + x®i®2] dt + v /tT dwi(t) , 
da\ = [—71a* + x^ia^] dt + ^/tT dwl(t) , 



da 2 = 
da% = 



X 2 , 
-72« 2 ~ —Q-x -\ 

'72«2 " 2 



dt-i 
dt 



72 dw 2 (t) , 
l2dw 2 {t) . 



(3.7) 



Here duik (t) is now a complex Gaussian white noise whose mean and variance are given by 

(dw k {t)) = , 
(dwk(t)dw*(t)) = Skidt . 



(3.8) 



The above equation is identical to the equation derived in positive P-representation when one discards the noise 
terms. This corresponds to the nonlinear classical equation for the OPO system. The main difference between the 
two sets of equations is the noise terms. In the semiclassical theory the noise is universal for all modes and comes 
from to the vacuum fluctuations, while in the positive-P equation the pump has a noiseless amplitude and the signal 
noise comes from the nonlinear coupling. However, we note that the Wigner equation after truncation is no longer 
completely equivalent to quantum mechanics, since it always leads to a positive Wigner function - thus, not all 
quantum states can be represented. 
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C. Observable moments and spectra 



The positive-P stochastic method directly reproduces the normally ordered correlations and moments, while the 
Wigner representation reproduces the symmetrically ordered moments. We also have to distinguish the internal and 
external operator moments, since measurements are normally performed on output fields that are external to the 
cavity. The technique for treating external field spectra was introduced by Yurke || , and by Collett and Gardiner M . 

These external field measurements are obtained from the input-output relations of: 



$f = . 2~!™%{t) - <J>f (t) , (3.9) 



j \-j - y - lj ™3 v v i 

where and 3>° ut (t) are the input and output photon fields respectively, evaluated at the output-coupling mirror. 

The most efficient transport of squeezing is obtained if we assume that all the signal losses occur through the output 
coupler, so that 71 = 7°"*. We will assume this to be the case. 

The crucial quadrature variables of the system have the definitions: 

Xj = (aj +atj , 

y j = ~(a j -af) . (3.10) 



1 

There are also corresponding external quadrature field variables, defined as: 

Xj = + $°" tf 

1 

i 

Similarly, we can define c- number stochastic quadrature variables within the relevant representations, thus giving: 



Xj = (otj +a+) , 

Vj = \ {<*j - at) (3.12) 

Of especial interest is Y\ since this is the low-noise, squeezed quadrature. Here we note that the instantaneous 
correlation functions of the intra-cavity field operators are called the moments. Typically, they are not easily measur- 
able, when compared to output moments or spectra, but they are useful in that they provide a check on the accuracy 
of the calculation of measurable spectra. 

The squeezing in terms of the intra-cavity quadrature variances corresponds to an instantaneous measurement of 
the field moments. If such a measurements were possible, it would include contributions from all frequencies. For 
measurements averaged over a long time T, it is the low frequency part of the spectrum that is the relevant quantity, 
and we shall focus on this, as it usually determines the maximum squeezing possible. The output measured spectral 
variance V® of a general quadrature 

can be written 

Vf(u)6(w W) = (AXfMAXjV)) . (3.13) 
where the fluctuations AJ* are defined as AXj* = Xj—(Xj) and the frequency argument denotes a Fourier transform: 

J V 27T 

Since the P-representation is normally ordered, it automatically provides the normally-ordered moments: 

(:x e j (t)x e j (t):)={x e j (t)x e j (t)) P . (3.14) 
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Also, the +P spectral correlations correspond to the normally ordered, time-ordered operator correlations of the 
measured fields. We therefore define Fourier components of the normalized quadratures as: 

x* (t) = / ^§e-^ H , (3-15) 

This leads to the following well-known result for the general squeezing spectrum, as measured in an external 
homodyne detection scheme: 

V?(u)5(l> +u') = 1 + 2 7 f * (Ax* (w) Axj (c/)) p • (3.16) 

Note that vacuum (input) field terms do not contribute directly to this spectrum, as they have a vanishing normally- 
ordered spectrum, and are not correlated with the coherent amplitudes in the +P representation. 

In the case of the Wigner representation, the correlations and moments are given with symmetric ordering. Thus, 

for example, (a*j (t) ctj (t)) w = Oj(t),aJ(t) /2^ = 1/2 in the vacuum state. The normally-ordered internal 

field moments are easily calculated, by using equal-time commutators to change the ordering from symmetric to 
normally-ordered : 

(:%(t)& j (t):) = (x 9 j (t) a *(t)) w -l. (3.17) 
Similarly, the normally-ordered squeezing spectrum, as measured in an external homodyne detection scheme is: 

Vf(Lu)S(uj + uj') = (X? M X\ (a/)) . (3.18) 

It is essential here to include the vacuum field contributions from reflected input fields, as these are correlated with 
the internal Wigner amplitudes, and hence have a significant contribution to the spectrum. In fact, these input fields 
can be shown to correspond directly to the noise terms in the relevant Wigner equations, leading to the identification: 

^ = V2$r(t) , (3.19) 

where $*™(t) is a c- number amplitude corresponding (in the Wigner representation) to the quantum vacuum input 
field. 

The fundamental property of the Wigner function is that the ensemble average of any polynomial of the random 
variable a and a* weighted by the Wigner density exactly corresponds to the Hilbert-space expectation of the cor- 
responding symmetrized product of the annihilation and creation operators. Therefore, the truncated theory with a 
positive Wigner function can be viewed as equivalent to a hidden variable theory, since one can obtain quadrature 
fluctuation predictions by following an essentially classical prescription; in which even the noise terms have a classical 
interpretation as corresponding a form of zero-point fluctuation. This cannot be equivalent to quantum mechanics in 
general, but may provide similar results to quantum mechanics under some circumstances. 

IV. BELOW-THRESHOLD PERTURBATION THEORY 

Next we wish to rescale the equations. This has the merit of showing explicitly how a small noise expansion can 
permit us to use a type of perturbation theory whose zero-th order solution is the classical solution, rather than the 
Feynman approach where the zero-th order solution is the free-particle case. In order to show this systematically, a 
formal perturbation expansion in powers of g is now introduced, where the scaling parameter g is given by: 



.g = l/ v /47e7i"=l/ v / 2iV c 7r , (4.1) 

where N c = 2I c /^2 is the threshold pump photon number, and a dimcnsionless decay ratio, j r = 72/71 is introduced. 
An equivalent definition is: 

V (4.2) 



V 2 7i72 



This clearly determines the ratio of nonlinear to linear rates of change. Next, we introduce a scaled time r = 71 1 
and a dimensionless driving field /i = £/£ c = I {ill?), so that the equations can be expressed in terms of the three 
dimcnsionless parameters g,n,j r . Finally, we expand the scaled coordinates in a power series in g , to give: 
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(4.3) 



The expansion given here has the property that the zero-th order term corresponds to the large classical fields of 
order 1/g, while the first order term corresponds to the quantum fluctuations of order 1, and the higher order terms 
correspond to nonlinear corrections to the quantum fluctuations, of order g and greater. For a given fundamental 
decay rate 71, the expansion coefficient g 2 is inversely proportional to the input photon flux required to obtain the 
threshold condition. Thus, the smaller g 2 is, the larger the required input field. 

A. Matched power equations in positive P-representation 

Here we will first be interested in the analysis of the steady state moments. Subsequently we will calculate the 
spectral correlations of the solutions using Fourier transforms of the calculation done in the time domain. The 
equations for the quadrature variables in the positive P-representation are 



dx\ — 

dyi = 
dx 2 = 
dy 2 = 



X 



-722:2 - - (x\ - yl) + 2£ 



dt + 

dt — i 
dt 



X 
X 



\fx 2 + iy 2 dwi(t) + \lx\- iy\dw 2 (t) 



y/x2+ iy 2 dwi(t) - yx\- iy\dw 2 {t) 



-72J/2 - -xiyi 



dt 



(4.4) 



The stochastic equations are now solved by the technique of matching powers of g in the corresponding time- 
evolution equations. This technique can be analyzed diagrammatically, and so can be termed the 'stochastic diagram' 
method |E0(|. The zero-th order solution is: 



dx[ 0) = 

dyl 0) = 
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dr 
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(4.5) 



These equations are the classical nonlinear equations for the cavity, expressed in terms of the quadrature amplitudes 
of dimcnsionless scaled fields. The steady-state solution below threshold is well-known, and is given by: 



..((') 



yi 0) = ,f = 0; 4 0) = 2 M 



(4.6) 



With no loss of generality, we can set all odd orders of X2,V2, and all even orders of 2^ , to zero, since one 
can set these to zero initially, and these orders do not change in time. To first order, the equations are given by: 

dx^ = — (1 — [i) x^dr + y^2fldw x (T) , 

' ' ) j (4-7) 



dy[^ = - (1 + fi) yi'dr - i^/2jidwy(T) , 
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where, dw x ^(T) = (dwi(r) ± cIw2(t)) /\/2. These equations are the ones that are normally used to predict squeezing. 
They are non-classical, but correspond to a very simple form of linear, non-classical fluctuation which has a Gaussian 
quasi-probability distribution. In other words, if no higher-order terms existed, the result would be an ideal squeezed 
state in the sub-harmonic, together with an ideal coherent state in the pump. 

Of more interest to the present paper, is the behaviour to the next order. This is the first order where nonlinear 
corrections to ideal squeezed-state behaviour will occur. We find the following: 



dxf = 



,.(2) 



- X 

2 



dr 



(2) i (1) (1) j 



(4.8) 



While we do not wish to include any effects beyond the first nonlinear corrections, it is not possible to consistently 
neglect the third-order in perturbation theory. This is because the first non-trivial correlations arise in terms like 



which have the same formal order as terms of the type ( 



x (3) x (i) 



) . Therefore, to obtain a consistent 



expansion for the correlations that are of interest, we must compute the third-order terms as well. These satisfy the 
following equations: 



(3) _ 



dx{ 



dyf ] = 



(1-^ + 5(^+^1 



^ __ L 4 2) du; x (T) + iy ( 2 z 'dw y (T) 
y% dw x (r) — ix% dw y (r) 



(2), 



dr 



dr 



2^ 



(4.9) 



The equations of this order have a non-trivial noise term, which depends on the second order pump quadrature 
solution. 



1. Operator moments 



We now wish to calculate the operator moments. To proceed further, we use ltd calculus to derive stochastic 
equations for quantities of interest, which in the present calculation are y^yi and y[^y^ . These equations contain 
quantities involving variables lower down in the hierarchy, as well as terms generated from the noise correlations. 
Finally, we compute the the steady state averages of the quantities of interest, so that the noise terms vanish. In the 
present case, this yields, 
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x (1) x (1) 
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7r 



( 7r + 2) \1 — fj, 



(4.10) 



The first quantity above is related to the depletion of the pump that supplies the energy for the sub-harmonic 
mode. The follow two quantities are the squeezed and enhanced quadratures normally obtained in the linearized 
theory, while the fourth one is the first correction to the linearized calculation. The last one is the steady state triple 
quadrature correlation. This quantity has been suggested previously as a way to test quantum mechanics against a 
local hidden variable theory pl| . 

The steady state intra-cavity squeezed quadrature fluctuations are obtained as: 
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(4.11) 



The intra-cavity squeezing quadrature near threshold is not perfectly squeezed, as the nonlinear correction is 
divergent near this point. This is shown in Fig. ([!]). It is clear that the nonlinear corrections to the overall moment 
scale as g 2 / (1 — fi) , and hence only give large corrections extremely close to threshold, with /x « 1 — g 2 . 




FIG. 1. Squeezing moment (j/i) versus driving field ft, with g 2 — 0.001, 7,. = 0.5 
Considerations related to optimal squeezing will be treated later, in the frequency domain. 



B. Matched power equations in semiclassical theory 

We can scale the quadratures variables in in semiclassical theory in the same way as before. Firstly, the equations for 
the quadratures are: 



-I1X1 + -z {xix 2 + 2/12/2) dt + sTU [dwi(t) + dw 1 (t)\ , 
X 



dyi = -712/1 + g (xiV2 - x 2 yi) dt-iy/y[[dwi(t) + dwl(i)} 
(x\ - yl) + 2£ 
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dx 2 = -72^2 - 7 (^1 - 2/i) +2£ dt + ^[yi [dw 2 (t) + dw 2 (t)] , 
dy 2 = ^-722/2 - 2 Xiyi \ dt ~ i ^ 2 \. dw ^\ t ) ~ rfw 2 (*)] • 
In the new scaled time, the correlation function of the noise terms is: 

<&(*)£*(*')} = (6 (t/ 7 i)C (r'/li)) = 7i«r r') = 7i(&(t)£?(t')> 



(4.12) 



(4.13) 







where we have written the Wiener increment as dw(t) = £(t)dt. Next, we redefine the white noise that drives the 
stochastic semiclassical equations as 



dw xl ( 2 )(T) 



dw 1{2) {r) +dw* 1{2) (i 



V2 
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(4.14) 



The dimensi onle ss driving field [i is introduced as before, and the Wiener increments dwi(r) have the same properties 
as defined in (3.8), except for changing t to the dimensionless scaled time r. Next, we use the same technique of 
matching the powers of g in the corresponding time-evolution equations. The zero-th order equations are: 
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y\ 2/2 



0),,(0) (0) (0) 



(0) (0) 



dr , 
dr , 
- 2/i 



dr 



(o) 

y 2 



T (0) ? .(0) 



dr . 



As in the positive-P case, the steady state solution below threshold is given by: 

J0) 



(o) (0) 



2/f = 



(4.15) 
(4.16) 



'-2 — 2/* ■ 

To first order, the equations are given by: 

dx^ = — (1 — /Li) x^dr + \/2dw x i(T) , 
dy[ 1] = - (1 + m) Vi ] dT + V2dw yl (T) , 
dx 2 = —"frX^ dT + 2^/ r dw x2 (T) , 
dy^ = -jryi^dr + 2-f r dw y2 (T) . 

While the zero-th order equations are essentially classical, in this first order set the noise appears as a quantum 
effect. This is still a linear approximation, as all nonlinear corrections come from the next orders. 
The second order equations are: 



(4.17) 



dy\ 



(2) 



(2) 



dx 



dy ( 2 2) 



(l-M)4 2) + 5(4 1) 4 1) +2/f^ 1) ) 



dr 
dr 



(2) 



X 



{1) y[ 1] 



dr . 



(4.18) 



We need to go beyond this order in perturbation theory to compute the first nonlinear corrections. The third order 
equations are: 



dx? } = 



A ( 3 ) 

dy\ = 



dx 2 = —1 
dy ( 2 3) 



- (i - /*) 4 3) + \ (4 X) 4 2) + 4 2) 4 X) + W + itfW) 

-(i + /*)yf ) + 5(*i 1) % ca) + «i : 



^U 1 ) T (!) 7 ,( 2 ) 

i/2 x 2 s/i 



(2) (1) 

4 2/i 



dr 
dr 



(1) (2) 

2/1 2/i 



(3) , (1) (2) . (2) (1) 

2/2 +^1 V\ + x \ 2/i 



(4.19) 
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1. Operator moments 



The steady state averages of the quantities of interest can now be using the truncated Wigner distribution, therefore 
obtaining the symmetrically ordered correlation functions: 



•r 



(2)\ 
2 / 



<vfV w > 

I (1) (3)\ 



/' 



l-/x 2 
1 



1 



7r 



7r 



2 V7r + 2/ (1 - fi 2 ) 2(1 + /i) 2 [ 7r + 2(1 + (Lt)] ' 

-7r , 7r(2-Ai) + 2(l + A*) ' 



4(l+ M )(l-/i2) 

7r 



7,. + 2 (l + /i)[ 7r + 2(l + /i)] 



1 



7 r + 2 / \ 1 - /i 
2 7r 



7,. + 2 y v 1 - m : 



1 



(4.20) 



The main difference in these calculation compared with the positive-P result, appears in the nonlinear correction 
for the sub-harmonic squeezed quadrature. Up to second order in g we have 



(vi) 



1 



2/ (1) (1)\ , 4/ (2) (2)\ . r, 4/ (1) (3)\ 



2(1 + M )(1- M 2 ) 



7 r 7 r (l + 2^-2^ 2 ) + 2/x(l + /i) 



7r 



(l + /x)[7r+2(l + /i)] 



(4.21) 



The similarities and disagreement between this result and the positive-P expression for the same quantity deserve 
further comments given in the conclusion section. In particular, we notice that while the linear term agrees, the 
nonlinear term is not in agreement well below threshold. 

This comparison is shown in Fig(0), which compares the nonlinear parts of the moment in the two representations. 
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FIG. 2. Nonlinear correction to the squeezing moment {^■Vi') versus driving field /x, with g 2 = 0.001, j r = 0.1, 1, 10 . Solid 
line corresponds to the positive-P representation, dashed line to the Wigner representation. Best squeezing occurs with smallest 
value of j r . 

Just below threshold both theories give nonlinear corrections which are essentially identical. There is also good 
agreement in the limit of -f r — > 0, where 72 <C 71; but for -f r > and driving fields below threshold, there is 
substantial disagreement in the nonlinear corrections to the squeezing between the two representations. This can be 
attributed to the neglect of third-order quantum correlations in the truncated Wigner representation, which results 
in the appearance of nonlinear squeezing effects even in the limit of zero driving field. Such effects are due to the 
scmiclassical vacuum inputs, which do not appear in the positive P-representation. 



V. SPECTRAL CORRELATIONS 



Next, we proceed to analyze the problem in the frequency space by taking the Fourier decomposition of the fields 
in order to understand the role of the first nonlinear correction in the squeezing spectrum. It is important to stress 
that most of the measurements performed are done in Fourier space. 

The nonlinear corrections to the spectrum have a strikingly different behaviour to the case of the squeezing moments. 
The reason for this is that the nonlinear corrections are due to low-frequency, narrow-band critical fluctuations. These 
have a very small effect on the moments, which correspond to an integral of the spectrum over all frequencies, unless 
extremely close to threshold. However, they can have a very large and disruptive effect on the very important zero 
frequency component of the squeezing spectrum, where the quantum noise is at its lowest level. 



A. Positive P-representation 

The spectrum can be calculated directly from the Fourier transform of the stochastic equations. We also represent 
the white noise that drives the stochastic equations by its Fourier transform £ XtV (O), where the spectral moments of 
the stochastic processes are: 

(& («)) - , 

(&(fi)&(n')) = s ab s(n + n') . (5.1) 
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It is also useful to introduce a standard convolution notation, where: 

[A*B\(Sl) = f ^=A(SV)B(Sl - O') . 
J y2n 

The stochastic equations may now be rewritten in the frequency domain as - 
• First order: 

-(l) , Q x _ VZftx (n) 
Xl (U) ~ (ifi + l-p) ' 

(1 ) _ iy/^y (O) 

2/1 - PTTTm) • 



(5.2) 



(5.3) 



• Second order: 



4 2) («) = - 



.(i) ,(i) ,(i) ,(i) 



7r 



-(1) 

^1 



2 (in + 7 r ) 
(fi) 



(if) + 7 r 



(5.4) 



Third order: 



5p } (fi) = 

y{ 3) (fi) = 



-(2) 
2/2 * 



2 (iO + 1 - ju) 



2 (iO + 1 + /i) 



(5.5) 



-Z. Squeezing correlation spectrum 
We now calculate the spectrum of the squeezed field yi, which is given by (yi (fii) yi (^ 2 )). Thus, we obtain 
(m (Ox) ft (fi 2 )) = (yP (00 yf } (Q 2 )) + 

+ y 2 (yf > (fi 2 ) y[ 3) (Oi) + [Oi <- fi 2 ]\ + • • • (5.6) 



The contribution from the first order perturbation theory is the usual linearized squeezing result, given in this case 
by: 



yf } (QJtf' (fi 2 ) 
Similarly, the complementary (unsqueezed) spectrum is: 



2^(5(^1 +fi 2 ) 



ft 2 + (1 + m)' 



2^(fii+ft 2 ) 



(5.7) 



^ 2 + (1 - m) 5 



(5i 



Also, we can obtain the next order contribution to the squeezing, by calculating (y^ (fii)y^ (^2)) • To check the 
results, we can compare with the moment calculations, since: 



(5.9) 
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Using these results, we find that the internal spectrum of the squeezed quadrature, to this order, is given by 



and the squeezing spectrum is calculated to be: 



fi 2 + (1 + ^) 2 [ft 2 + (1 + m) 2 ] 2 

(o 2 + i- m 2 ) (i - M + 7r )(i + /i) -n 2 



2/i 7r (i-/i 2 ) (i-n)[n 2 + (i-v + lr ) 2 ] 

(1 +/i + 7r)(l + M) - ^ 2 



(l + /i)[fi 2 + (l + A i + 7r ) 2 

The corresponding external squeezing spectrum is then: 



V(to) = 1 - 



2 2 
M 7r 



n 2 + (i + /i) 2 [r> 2 + (i + Ai ) 2 ] 2 
(o 2 + i-/i 2 ) , (i- M + 7r )(i + M ) 



n 2 



2/i 7r (l-/i 2 ) (l- M )[n2 + (l- M + 7r ) 2 ] 

(l+/X + 7r )(l + /i) - fl 2 1 



(l + ^)[ft 2 + (l + ^ + 7r ) 2 



(5.10) 



(5.11) 



(5.12) 



This equation gives the complete linear and nonlinear squeezing spectrum, including all the nonlinear correction terms 
that contribute to order g 2 or 1/N. An illustration of the behaviour of the total spectrum is given in Fig 
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FIG. 3. Total OPO squeezing spectrum with g 2 = 0.001 , 7r = 0.5. The /i values plotted are /i = 0.1,0.3,0.5,0.7,0.9; larger 
values of \x give the most squeezing (lowest spectral variance). 

Fig (jij) shows how the nonlinear contribution changes with driving field, giving just the portion of the spectrum 
proportional to g 2 . 
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FIG. 4. Nonlinear OPO squeezing spectrum with g — 0.001 ,j r = 0.5. The maximum fi value plotted is fi — 0.95 



2. Triple Spectral Correlations 

Next, we can calculate the triple spectral correlations, giving as in the moment calculations: 
(ii(fii)jfi(fi 2 )»2(n3)) = g(x^(n 1 )y[ 1) (Q 2 )y 2 2) (fi 3 )) . 

Solving for , we nave 



,-,(2) 



(ifi 3 + 7 r ) 

Substituting from the first order spectrum, the final result to this order is obtained to be: 





fi 2 + n 3 ) 


(zf7 3 + 7 r ) 


^ + (i-m) 2 " 




o| + (i + M ) 2 " 



(i« (fiOjf (fi 2 )y^ (fi 3 )) = 



To check this result, we can evaluate moments: 
On integrating, we obtain the same result as in our moment calculation, given above. 



(5.13) 



(5.14) 



(5.15) 



(5.16) 



B. Semiclassical theory 

We will now compare these results with the corresponding results calculated in the semiclassical theory. Some 
differences between them could be an interesting test comparing quantum mechanical predictions with a hidden 
variable theory. 
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Again, the spectral correlations are calculated from the Fourier transform of the stochastic equations. In the 
frequency domain, the equations are written as 



• First order 



• Second order 



4 1} (n) 



v^gxi(n) 

V2g„i(n) 
(in + i + n) ' 

27r^x2(») 
(id + -f r ) 

27r^2(») 

(in + 7 r ) 



(5.17) 



s£ 2) (fi) 
yf(fi) 



4 2) (o) 



yf (fi) 



-(l) -fi) . -(l) -fi) 



2 (ifi + 1 - n) 



'2/2 



2/i * »2 



7r 



2 (ifi + 1 + /i) 

.(i) ,(i) .(i) ,(i) 
x\> * ^ ' - y} * y\> 



(fi) 



2 (in + 7 r ) 



(fi) 



(in + 7 r ) 



(5.18) 



Third order (sub-harmonic field) 



4 3) (n) 



yf>(n) 



-(1) -(2) . -(2) -(1) . -(1) -(2) . -(2) -(1) 

x\ * x\> + x\> * 4 + yj ' * y 2 ; + y| * y\> 



(fi) 



2 (in + 1 - n) 

x^ * y^ + * y^ — yj 2 ^ * a;^ — yj 1 ' * a:^ 



(fi) 



2 (in + 1 + /i) 



(5.19) 



i. Squeezing Correlation spectrum 
The spectrum of the fields are given, for instance for the squeezed quadrature j/i, by 

(Si(fii)Si(n 2 )) = (^(noy^^)) + . 9 2 {(y( 2) (n 1 )yf ) (n 2 )> 

+ (y^nOyf } (n 2 )) + (^(O^f^O)} + • • • (5.20) 
The first order perturbation theory generates the usual linearized squeezed result as in quantum theory 

and, similarly, for the amplified fluctuation quadrature 

<»i 1) (ni)4 1) (" a )) = ^; i + _^ . ( 5 - 22 ) 



16 



and for the pump quadratures, there is no first order squeezing 



(x£\sii)xP{Sh)) = (vPWvPiSh)) = T^P^s(n 1 + n 2 ) . (5.23) 

a 1 -f- 7 r 



The next contribution to the squeezing field quadrature is: 



(vf ) (n l) #(o 2 ))= >/^ +0 ^ 



1 - fl + 7 r 



n\ + (1 + M )2 1 (1 - M ) [n? + (1 - M + 7 r) 2 

1 + /i + 7 r 

(1 + //) [n 2 + (1 + /i + 7r )2 



and 



x < - 



[n; + (i + /i)s 
(i + /i)(i-/i + 7r )~nf 
(1 - M ) pi? + (1 - n + lr ) 2 ] 
(i + /i)(n-M+7r)-n? (i + m) 



(l + M )[n2 + (l + /i + 7r )2] 7r (l- M 2) 

and the internal (symmetrically ordered) squeezing spectrum is 

Q(O) = 2 + ffV f Ml+M) 

^ ' ^ 2 + (l+/i) 2 + [n2 + (l+^)2] 2 l7r(l-/^ 2 ) 

(1 - M + 7r)^ 2 + [(1 + M) 2 + 2yu(l + m)] (1 + M + 7r) 
+ (T+m) P 2 + (l + M + 7r) 2 ] 

(1 + ^ + 7r )fl 2 + (1 - /i 2 )(l - /X + 7r ) 
(l-/i)[fi 2 + (l- M + 7r ) 2 ] 



(5.24) 



(5.25) 



(5.26) 



Of greater interest is the external squeezing spectrum, which is obtained by including both internal fields and the 
correlated reflected vacuum noise terms: 

V(Q) = 1 - 4 M 2g 2 7r f M (l + r> 2 - M 2 ) 

^ ^ tt 2 + (l+/i) 2 [fi2 + (l +M )2] 2 \ 7,(1 "/**) 

[(1 - /i)(l - /i + 7 r ) - 2/i 2 ] ft 2 + [1 - fl + J r ] [1 + fi + V 2 + ^] 



(1-/*) [fi 2 + (l-M + 7r) 2 ] 

[(l + ^)(l + /i + 7 r ) + 2/i 2 ] 2 + [l + ^ + 7r ] [l + 3/i + m 2 - M 3 ] 



(l + /i) [fi 2 + (l + ^ + 7r) 2 



(5.27) 



This semi-classical spectrum is quite different from positive-P calculation when fi — ► but gives a compatible result 
near threshold, that is in the limit fx — > 1. A detailed comparison of the zero- frequency behaviour is shown in Fig (|^). 
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This means that even when the pump is off, semiclassical theory gives a distorted vacuum spectrum due to the 
presence of the nonlinear crystal. This happens because in this theory the vacuum fluctuations are taken as real, 
and then two vacuum modes can interact inside the crystal as real fields. In the limit of 7r — > , the two spectra 
become compatible again, as the semiclassical theory decouples the second-harmonic mode from its vacuum input in 
this limit. In the case of threshold fluctuations, we can interpret the agreement as due to the large photon numbers 
involved - which means that the truncation approximation used for the semiclassical calculation is more reliable. 



C. Optimal Squeezing 

It is interesting to evaluate the squeezing or low-noise quantum correlations in the limit of zero frequency, that is in 
the resonance regime which is generally the frequency of maximum squeezing. We obtain from the positive-P result: 



V(0) 



4 M 



2^ 2 



1 



4 7rM 2 ( 7r + 2) 
(1 -/*)[(! + 7r) 2 - A* 2 ] 



(5.28) 



Near threshold, where fx w 1, we can set fx = 1 + 5, and expand in powers of 5 < 0. Minimizing this result with 
respect to 5, we find that, to leading order in g, the optimal driving field is the solution to the following equation: 

S 3 (2S + 7r ( 7r + 2)) 2 = 5 2 7r ( 7r + 2)(4J - >( 7r + 2)) 

This is a quintic equation, but it has simple closed form solutions in two limits, depending on whether 7r 3> g 2 / 3 
or 7r <C g 1 ! 3 . In the first case, the variance can be rewritten as: 



V(0) = 



1 



2.9 2 



(5.29) 



this result with respect to 5 , we find that the minimum level of internal fluctuations occurs in a narrow 
frequency range near ft = , at a driving field just below threshold, with S = —g 1 ! 3 so that: 



f^opt 



.2/3 



(5.30) 
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To leading order in g , the corresponding spectral variance is: 

v oP t(o) = - A g i/3 



(5.31) 



This result of an N~ 2 / 3 scaling confirms an approximate calculation of Plimak and Walls JTJ], although the self- 
consistent method used by these authors makes it difficult to obtain the relevant driving field. 

The physics of this is clearly that the onset of critical fluctuations starts to spoil the noise-reduction even before the 
critical point is reached at fj, = 1. For example, with 7 r « 1 and /i — 0.9 , we find that V(0) ~ 0.7 x 10~ 2 , or about 
2\dB below shot noise, as predicted from the analytic theory. This can also be seen from the way that the third order 
term includes contributions from the critical fluctuations in x\. A direct calculation from the full spectrum shows 
that this is a true minimum for all frequencies, even including Q > . 

However, the situation clearly changes as 7 r — > , in which case much greater levels of spectral noise-reduction are 
possible. This is plotted in Fig M: 
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o 
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0.95 



1 
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FIG. 6. Optimization of zero-frequency squeezing spectrum versus driving field using the positive-P method, with g 2 = .001. 
Values of 7,. = 0.001, .01, .1, 1., 10. are used for the different lines plotted, with the lowest values of y r giving the best results 
for squeezing. 



Analytically, this limit gives the following result, provided that g 2 <C j r <C g 



2/3 . 



V(0) = \ 



2 I 2g*jr 
2 S 2 



(5.32) 



Minimizing this result with respect to S , we find that the minimum level of internal fluctuations occurs at a driving 
field very close to threshold, with: 



Mo pt = l-.9 1/2 (27,) 1/4 



(5.33) 



The corresponding variance is therefore: 



P opt (0) = . 9v ^72«5 4/3 ■ 



(5.34) 



This result can be much smaller than predicted by the calculation of Plimak and Walls [Q, since the damping 
ratio can be reduced (at least in principle) to an arbitrarily low level - although still bounded below by g 2 , in order 
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for perturbation theory to be applicable, so that we do not expect to obtain V op t(0) < g 2 . Of course, there are 
experimental limitations on this, due to absorption losses in the nonlinear medium at short wavelengths. Thus, for 
example, with the same value of g 2 = 0.001 as previously, but with j r = 0.01 , we find that the minimum spectral 
noise is predicted to occur at a driving field of /i = .93 , with a squeezing variance of 2.2 x 10~ 3 , or about 27dB 
below shot noise - about 6dB lower than before. 
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FIG. 7. Optimal zero-frequency squeezing spectrum versus frequency using the positive-P method, with g 2 = .001, ~f r = 0.01. 
Driving fields of /i = 0.9, 0.93, 0.96 are used for the different lines plotted, with the higher driving fields giving the best results 
for squeezing, except at zero frequency. 

This operating regime also has the property that the optimum frequency of noise reduction moves away from zero 
frequency as the driving field is increased above the optimum value, towards threshold. At slightly higher driving 
fields than the optimum point, a bifurcation to a spectrum with two minima occurs, although with similar levels of 
noise reduction, as shown in Fig (|^) . In this regime the results of the perturbation theory need to be checked by a full 
simulation of the stochastic equations. We have carried this out (see next section), and find that the full simulations 
do agree very well with the analytic predictions, even with this small damping ratio. 

D. Numerical Simulations 

The value of the nonlinear correction to the spectrum of the squeezed quadrature V(f2) can be worked out from 
a full numerical simulation j2j| of the relevant nonlinear stochastic equations. The optimal squeezing in the zero 
frequency part of the squeezing spectrum is predicted to scale as ./V -2 / 3 with roughly equal values of decay rates. For 
the simulations, we chose values of TV = g~ 2 = 10 3 , 7 r = 0.5. The simulations used a total dimensionless time-interval 
of T max — 1000. To ensure equilibrium, only the last 500 time units were utilized in the Fourier transforms. Time steps 
of At = 0.1 and At = 0.2 were compared to ensure convergence. The algorithmic technique is described elsewhere 
, and uses a semi-implicit central partial difference technique. To obtain the small nonlinear corrections near the 
optimum squeezing, we simulated the difference between the linear and nonlinear forms of the stochastic equation, in 
order to minimize sampling errors. It was also useful to initialize the x quadratures with a Gaussian ensemble close 
to the known steady-state variance, in order to reduce the time taken to achieve equilibrium. Typically, the relative 
error in the correlations due to finite step-size was around 10 -4 with these step-sizes. 

For these parameters the optimal driving field is predicted to occur at /i = 0.9, or approximately 80% of the critical 
intensity. We used 10 5 trajectories to improve the relative error due to sampling with a finite trajectory population, 
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giving relative sampling errors of less than 10 
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FIG. 8. Numerically simulated optimum nonlinear squeezing with g 2 = .001, j r = 0.5, fi = 0.9. Solid line is the +P 
simulation result, dashed line the analytic prediction from perturbation theory, dashed-dotted line the Wigner prediction. 

The calculated squeezing moment from the SDE simulations was: (yf) + 0.5 = .0271 ± 10 -4 . This is in excellent 
agreement with the below-threshold expansion, which gives (Yj 2 ) + 0.5 = .0272, as this is just outside the critical 
region. 

We find that the spectral predictions are also well verified by the simulations. These resulted in a value for the 
nonlinear correction to the zero-frequency spectrum, of AV(0) = V(0) — V^(0) = 3.75 x 10~ 3 ± .02 x 10~ 3 . By 
comparison, the analytic theory, worked to fourth order in g, gives the prediction that AV(0) — 4.02 x 10~ 3 . The 
residual difference of about 5% - which is significant compared to sampling error - can be attributed to the fact 
that there are higher order corrections that are not included in the analytic theory, and these are more significant 
in the zero-frequency spectrum than they are in the moment calculation. Fig. (||) shows the detailed results of the 
simulation. 
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FIG. 9. Numerically simulated optimum nonlinear squeezing with g 2 — .001, 7,. = 0.01, jj, = 0.93. Solid line is the +P 
simulation result, dashed line the analytic prediction from +P perturbation theory, ( results are identical with the Wigner 
prediction). Analytic predictions in this case are very close to numerical simulation results. 

In the analytic theory, we found that a smaller decay rate for the second-harmonic is predicted to yield a better 
squeezing optimum, as a function of driving field. In Fig. (Q) we verify this to be the case, by carrying out a full 
numerical simulation for j r = 0.01, /1 = 0.93, with T max = 2000, time steps of At = 0.05 and At = 0.1 for error- 
checking, and 10 4 trajectories. The results show that the analytic predictions and numerical simulations are almost 
indistinguishable in this regime. The sampling error was relatively larger, possibly due to the fact that the absolute 
noise levels are lower here. The agreement indicates that the perturbation theory is an excellent approximation to 
the full nonlinear equations with these parameters. 

VI. CONCLUSION 

We have calculated the nonlinear quantum fluctuations in a parametric oscillator below the classical threshold, using 
a nonlinear stochastic positive-P theory, with both asymptotic approximations and a numerical technique. There is 
excellent agreement between numerical and analytic calculations. Corresponding results for the Keldysh diagram 
method require a summation over infinite sets of diagrams, in order to fully include the reservoirs. The advantage of 
the present method is due to the fact that the coherent state basis is a more natural basis set for an open system, 
since it allows the damping reservoirs to be treated non-perturbatively. 

Optimal squeezing in the output spectra corresponding to these moments were estimated. We found that the best 
squeezing in the zero frequency part of the squeezing spectrum scales like TV -2 / 3 just below threshold, provided the 
two fields have similar damping rates. In other words, at the true critical threshold - where the linear squeezing is 
optimized - the nonlinear corrections are too large to give the lowest overall zero-frequency squeezing. Instead, one 
should operate below the critical point to optimize the spectral squeezing. Using an entirely different method, a 
calculation by Plimak and Walls (l4| also predicted that the optimum zero frequency squeezing spectrum scales like 
N~ 2 / 3 , or equivalently, as /~ 2 / 3 for a given input flux / . Our general scaling results agree with theirs, except with 
a different spectrum. We attribute the difference to the systematic +P stochastic diagram procedure used here to 
calculate the spectrum, rather than the Feynman diagram method - which involves additional approximations. 

We also found a new regime in which the lower limit to the spectral noise-reduction depends on the decay rate of the 
second-harmonic field, which can be reduced to an arbitrarily low level. This has a reasonable physical interpretation, 
since the second-harmonic losses are essentially parasitic losses, which do not contribute to the desired squeezing 
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output. The ultimate limit to squeezing in this regime is set by even higher order terms in perturbation theory. 
We conjecture that optimization of both the driving field and the relative decay rate may result in a final squeezing 
variance scaling as Ji/I c ■ 

A calculation with the truncated Wigncr method, or scmiclassical technique, was also carried out. Well below 
threshold, we found that while the linear terms agreed with full quantum calculation, nonlinear corrections and 
higher order correlations tended to disagree, especially for high second-harmonic losses. However, near the critical 
point, the situation changed. Here, even though the dominant terms are nonlinear, we found excellent agreement 
between the two methods. 
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